Quasi-static crack propagation in heterogeneous media 
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The dynamics of a single crack moving through a heterogeneous medium is studied in the quasi- 
static approximation. Equations of motion for the crack front are formulated and the resulting 
scaling behaviour analyzed. In a model scalar system and for mode III (tearing) cracks, the crack 
surface is found to be self affine with a roughness exponent of f = 1/2. But in the usual exper- 
imental case of mode I (tensile) cracks, local mode preference causes the crack surface to be only 
logarithmically rough, quite unlike those seen in experiments. The effects of residual stresses are 
considered and found, potentially, to lead to increased crack surface roughness. But it appears 
likely that elastic wave propagation effects may be needed to explain the very rough crack surfaces 
observed experimentally. 



PACS numbers:62.20.Mk, 03.40.Dz, 46.30.Nz, 81 
Since the work of Mandelbrot, Passoja and Paullay [jj] , 
many experimental |^-|^] and numerical || studies on ge- 
ometrical properties of crack surfaces and crack fronts 
in heterogeneous materials have demonstrated their self- 
affine scaling properties. Cracks, in a large variety of 
materials from aluminum alloys Q to ceramics |3| and 
rock B , seem to leave behind surfaces whose height vari- 
ations are characterized by a roughness exponent £ ~ 0.8. 
In the meantime, there have been significant advances in 
the understanding of the dynamics of interfaces and lines 
pinned by quenched random impurities |^-|8| . Potentially, 
these might shed light on the problem of crack surfaces 
and fronts since a crack surface can be viewed as 
the trace left behind by the crack front as it traverses 
the sample. But such theoretical analyses of crack sur- 
faces have so far neglected the crucial long range effects 
of elasticity. In particular, the deformations of a crack 
front change stresses at distant points on the front flio| ]. 
Moreover, the stresses are also modified by the bound- 
ary conditions on the rough crack surface, which give rise 
to dependence on the history of the crack front. As the 
crack moves, it relieves both the elastic stresses generated 
by the external loading and those due to the presence of 
quenched impurities such as inclusions, micro-cracks or 
dislocations. Because of elasticity, these defects give rise 
to long-range correlated randomness in the equations of 
motion. 

All of these effects influence the macroscopic scaling 
behavior of the crack surface and have to be included in 
a realistic model of crack propagation. Unfortunately, 
the full elastodynamic problem is very difficult. Indeed, 
even the criteria that determine the local direction of 
advance of crack front are not well understood Jll|] . 

In this paper, we present the results of a study of a 
crack propagating through a heterogeneous medium in 
the quasi-static approximation. We first consider a sim- 
plified scalar elasticity theory and then the real case of 
vectorial elasticity. We use the symmetries of the sys- 
tem and elasticity theory to determine phenomenological 
equations of motion. 

Our conclusions are, unfortunately, that for the prin- 
cipal situation of experimental interest of tensile (mode 
I) cracks [and also shear (mode II) cracks] the predicted 
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roughness of the crack front is only logarithmic rather 
than a power law of the length scale. Tearing (mode III) 
cracks are substantially rougher; however, these tend not 
to be stable experimentally and thus the applicability 
of the results is questionable. In all cases, random long 
wavelength residual stresses can increase the roughness; 
the effects of these are discussed briefly. 




FIG. 1. Crack propagating through a heterogeneous 
medium. The x co-ordinate of the front is f(z,t) and the 
free crack surfaces S+,S- are located at y — h(x,z). The 
local tangent, forward and normal vectors t,a, n at a point 
vcf on the front are also shown as are the directions of mode 
I (tensile), mode II (shear) and mode III (tearing) loading. 

The geometry of the crack is shown in Fig. |l|. The 
crack front is oriented along the z— axis and moves in the 
positive x direction. The x (in-plane) coordinate of the 
front at time t, x c (z,t) — f(z,t) is assumed to be single 
valued in z. As the crack moves, it leaves behind free 
surfaces 6>+ and S- on either side of the crack, located in 
body coordinates at y = h(x, z) for x < f(z, t). Thus, the 
instantaneous position of the front is given by the curve 
tcf(z, t) = f(z, i)x + h[f(z, t), z]y + zz. We assume that 
the crack motion is quasi-static, so that at each instant, 
the scalar displacement field u(x, y, z) satisfies Laplace's 
equation with the boundary condition that the normal 
stresses on the surfaces 6>+ and S- vanish. The sample 
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is loaded far away from the crack, such that 



K°°y Ufa 



(1) 



for y/ ' x 2 + y 2 large, with K°° the applied stress inten- 
sity factor and M the elastic modulus. Displacement 
and stress fields in the medium are determined by the 
position of the front and the shape of the crack surface 
and these in turn influence the subsequent motion of the 
crack. We now derive general equations of motion for the 
time evolution of the crack and analyze them to deter- 
mine scaling exponents that characterize the roughness 
of crack fronts and surfaces. 

A local coordinate system can be defined at a point 
tcf on the crack front, with axes along the tangent t = 
d 7 rcF /\d z rcF\, forward direction a = d 2 rcF/\d 2 rcF\, 
and normal n = t x a (see Fig. |l|). The direction in 
which the crack front moves should be governed by the 
angular dependence of the stress and the local yield stress 
near the crack tip. The stress field close to the crack tip 
has the form 



K{r CF ) 



cos(0/2)n - sin(0/2)a] + S a (r CF )a + 0(y/p), 



(2) 



with p the distance to the crack front and 9 the angle 
with respect to a in the (a-n) plane. The local stress 
intensity factor is K(tcf) and S is the finite part of the 
stress whose normal component S n must vanish. The 
divergent part of the stress field does not break the re- 
flection symmetry with respect to the crack front plane 
(spanned by t and a), so the crack should move in the 
forward direction a leaving behind a crack surface that 
is smooth on the smallest length scales. 

As the crack front advances, an elastic energy 

G(^°°,rcF,{/},{/i}) = R2( m F) P er unit area of gen- 
erated crack surface is released by the medium. In or- 
der to create the new crack surface, a surface energy 
T(tcf) per unit area is required. This is a local material 
property. We will assume that the zone with non-linear 
deformations, plasticity or microscopic side branches is 
small and that linear elasticity theory is applicable out- 
side of this region; T will include the energy to break 
bonds, as well as that dissipated as a result of plastic 
flow etc. Under quasi-static conditions, the forward ve- 
locity v a = a • dtXcF of the crack front is proportional to 
the excess elastic energy released, i. e., 

^vaivcp) = G(K°°, v CF , {h}, {/}) - r(r CF ), (3) 

where \i is the effective mobility of the crack front. If fi 
is much smaller than (cT) , with c the speed of sound, 
the quasistatic approximation should be good except on 
very long length scales. 

Energy considerations alone cannot determine how the 
direction of motion of the crack front changes. Instead, 
we expect this to be determined by the small scale pro- 
cesses that break the reflection symmetry about the local 



crack surface. In particular, the finite part of the stress 
field near the crack front, S (1 (jccf) in Eq.(^) is antisym- 
metric about the crack surface. The competition of this 
with the local random variations in the material yield 
stress will determine the curvature n ee n- <9 a a of the path 
of the crack front as it advances: S a will tend to make the 
crack curve in the direction in which the "hoop" stress 
eg or the total stress |<r| is larger. For small changes in 
the direction of advance, this will be determined by the 
stresses near 9 = 0; from the small 9 behaviour of Eq.(|j), 
we thus see that k should be proportional to —S a . But 
the yield stress ( or other appropriate fracture stress) of 
the material, ay, can vary with position and the crack 
will tend to bend in a direction in which ay is lower. 
We postulate that the competition between these effects 
determines the local curvature: 



K = -(S a + bd n (Jy)/\, 



(4) 



with b~K 2 /ay 2 M 2 of order the distance from the crack 
tip at which the stress is ay and \^K c \/b, with K c the 
critical stress intensity factor for crack advance. The de- 
tailed form of Eq. (|4|) will turn out to be unimportant at 
long length scales. 

Both S a and the energy release rate, G, can be written 
as a sum of two terms. The first terms, G u and S% are 
the energy release rate and the nonsingular stress near 
the crack front for a homogeneous medium with the same 
loading and geometry of the crack surface, and thus do 
not depend on the heterogeneities in the system. These 
are functionals of the entire shape of the front and the 
crack surface, as well as the applied load. The second 
terms, G r and S^ are the contributions from the relief of 
frozen-in stresses due to the heterogeneities and defects 
in the medium as well as any elastic inhomogeneities. 

In general, G r and will also depend on the crack 
shape. But for almost straight cracks , they will primar- 
ily depend on the unperturbed shape of the crack and 
to linear order in the randomness this effect can be ne- 
glected. We can now define random variables 

X(x, z) ee bd n a Y [x, h(x, z),z] + S^[x, h(x, z),z], (5) 
j(x, z) ee r[cc, h(x, z), z] — r — G T [x, h(x, z), z], (6) 

on the crack surface, with (G r ) — (S^) — and (r) = T, 
the mean fracture toughness. These will depend on 
h(x, z) but to lowest order in h this effect is not impor- 
tant. The means of 7 and \ are defined to be zero and 
their correlations are taken to be 

(7(x,z)j(x',z'))=A(x-x',z-z% (7) 
(X(x,z)x(x',z')} = r(x-x',z-z'). (8) 

We first consider short range correlations where A and 
T oc 8{x — x')5(z — z'),with some short distance cut-off, 
later returning to the effects of long range correlations. 
The loads and the shape of the crack surface determine 
the stress field. Consequently Eqs.(|[]|) determine the 
evolution of the crack front. 
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We have obtained a solution for the displacement 
field perturbatively around a flat, straight crack in an 
isotropic, homogeneous medium that can be carried out 
order by order in h(x, z) and f(z,t) [12]. This solution 
combined with Eq.(|3j) yields an equation of motion for 
f(z,t). To linear order, out-of-plane deviations of the 
crack surface do not contribute to the divergent stress 
just ahead of the crack tip that determines the energy 
release rate. The behavior of the in-plane displacement 
/ is identical to that of a crack which is restricted to 
move in a plane: 



Hd t f{z,t) 



dz 



,f(z',t)-f(z,t) 



7T \Z — Z' 

- 7 (/(z,t),z) + G 00 -f, 



(9) 



with G°° = (K°°) 2 /M and V denotes the principal part 
of the integral. Similarly, Eq. (||) gives us an equation for 
the evolution of h{x,z). Stresses generated by in-plane 
fluctuations of the crack front do not break the reflection 
symmetry and therefore cannot affect the curvature k. 
Thus, the leading correction to S% involves only out-of- 
plane fluctuations h. The linearized equation for the evo- 
lution of h as the front passes a point tcf at x = f(z, t) 



oo oo 



Xd 2 h{x, z) = VG™V / dx' / dz'J(x -x',z- z')h(x' , z) 



(10) 



-oo — oo 



-X(x,z), 



independent of /. The Fourier transform of the long 
range kernel J is 



J{q,k) = {q-ieff 2 Y{k/ q ), 



(11) 



with Y a scaling function. To linear order the equations 
of motion for / and h decouple. Therefore the crack 
surface is determined independently of the dynamics of 
the front, which is not true in general. 

Equation has been studied both analytically and 
numerically 1 12| . The front is arrested for small external 
loads, and there is a critical load K^° (corresponding 
to G£°) at which the crack just begins to move. For a 
load K°° slightly above this threshold, the average ve- 
locity of the front scales as v ~ (K°° - K^Y . The 
motion of the front is rather jerky with fluctuations in 
the velocity correlated up to a distance £, which diverges 
at threshold like £ - {K°° — K™)~ 1/ . At length scales 
smaller than £, the in-plane front profile f(z,t) is self- 
affine with ([/(«, t) - f (z',t)] 2 )~\z - z'\ 2( t . A renormal- 
ization group e-expansion and numerical simulations sug- 
gest Q = 1/3, v = 3/2, (3 w 7/9. At length scales larger 
than £ (or well above threshold), the crack front moves 
more uniformly, and its in-plane fluctuations scale as 



([f(z,t)-f(z',t)} 2 )~log\z-z'\, 



(12) 



as can be seen by Fourier transforming Eq.(g). 

In order to calculate scaling properties of the crack 
surface, we next analyze Eq. (jToj) . At length scales larger 



than b, the left hand side of Eq. ( 10 ) becomes negligible 
and the direction in which the crack moves is determined 
by the competition between the non singular stress and 
the material properties near the crack front. Thus, be- 
yond this length scale our results will be the same if we 
alternatively assume that the crack tip moves in a di- 
rection in which the stress ag exceeds the local material 
yield stress furthest from the crack front. The crack sur- 
face roughness thus has the scaling form 

([h(x, z) - h(x', z')] 2 ) ~ \z - z'\ 2 ^H (\x - x'\/\z - z'\) , 

(13) 

with H a scaling function. The crack surface is 
anisotropic but with the same roughness exponent of 
C/j = 0.5 in both the x and the z directions. 

We now consider the real case of vectorial elasticity. 
Near the crack front, the stresses have the form 
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{K0 tj {9) + tfnSgO?) + K m J™(ff)}, (14) 



where the stress intensity factors K\,Kn and Km are 
associated, respectively, with discontinuities across the 
crack surface of displacements in the n (tensile loading) , 
a (shear loading) and t (tear loading) directions; and 
the £f„- are universal functions. If the crack is loaded far 
away purely in mode III, then out of plane z-independent 
deformations h{x) in the crack surface will not mix in the 
other modes near the crack front. The direction in which 
the crack progresses will then, as for the scalar case, be 
determined by the finite non-singular parts of the stress 
near the tip. The roughness exponents should thus be 
the same as the scalar case although the scaling functions 
will be different. Unfortunately, mode III cracks tend to 
be unstable E3], so the applicability of these results is 
questionable. 

The primary situation of experimental interest is mode 
I (tensile) loading. In this case, once the crack wanders 
out of plane, the local K\\ becomes non-zero, and with 
z-dependent distortions of the crack, so does Km. To 
linear order, K\\ and Km are functionals of h while K\ 
only depends on /. As the crack moves, the change in 

K 2 

the energy release rate G = ^l^jj-, where are the 
appropriate elastic constants Q, is dominated by 8K\. 
Therefore, as in the scalar case, the equation of motion of 
the in plane displacement will be independent of the out 
of plane displacement. The roughness of the front should 
be similar to the scalar case, but the out of plane rough- 
ness is very different. Following various authors, with the 
expectation that the local dynamics will be determined 
by the dominant terms consistent with symmetry, we as- 
sume that the crack tip locally prefers mode I loading 
|]l5| , ^6| . This implies that the direction of motion of the 
crack tip can be found by requiring that the curvature 
at the crack tip, k be proportional to the angular deriva- 
tive of the age component of the stress tensor, the hoop 
stress just ahead of the crack tip. As in the scalar case, 
this just makes sure that the crack is smooth on small 
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length scales, and the results on large length scales should 
be independent of the specific form of this condition. To 
linear order this leads to the equation for the out of plane 
displacement 



Vb 



Vb 

[Kf°d x h + K£\- X , 



(15) 



where K^, which is linear in h 



is the local mode II load- 
ing in the original x,y, z coordinates, obtained from the 
perturbative solution for the displacement field. Again, 
b is a microscopic length, \~K£^/b, and x is a random 
variable that includes the effects of local variations in the 
yield stress and the relief of random residual stresses. 

Once again / appears only as an implicit variable and h 
is determined independently of /. The Fourier transform 
of Eq. (|l5| ) involves a long range kernel 



J I (q,k) = \q\Y I (k/q). 



(16) 



We now see a crucial difference from the scalar case: J\ 
differs from J by a factor of l/lql 1 / 2 , thus the out of 
plane "stiffness" is much stronger in the mode I case due 
to the local mode selection. The extra stiffness leads 
to only logarithmically rough crack surfaces, in striking 
contrast to what is observed in experiments. 

A study of the non-linear terms that we have neglected 
throughout, leads to the conclusion that they will be ir- 
relevant at long length scales for all the cases considered. 
In our quasi-static approximation it thus appears that the 
only way to get rougher surfaces is if the random prop- 
erties of the solid have long range correlations, which we 
discuss next. 

Residual stresses cr r (r), which are present in the ma- 
terial before it cracks, will change the local stress inten- 
sity factors at the crack front for a given crack geome- 
try. If ((T r (r)cr r (r'))^l/|r — r'| Q at long distances, dimen- 
sional analysis yields additional correlations in the re- 
sulting G r and S£, ofthe form A(x,z) ~ |z| 1_a A(|iE/«|) 
and T(x,z) ~ |z| _a Y(|x/z|). Randomly distributed fi- 
nite sized defects such as heterogeneities, micro-cracks 
or loop dislocations each yield stress fields falling off as 
1/r 3 [[l7| and hence a — 3. For the roughness ofthe mov- 
ing front on scales larger than f , and the roughness ofthe 
crack surfaces in mode I , this is a marginal perturbation 
which will result in changing the logarithmic roughness 
to (In) 2 . Near the threshold for the crack motion, such 
correlations are also marginal. But in some disordered 
materials, growth processes may result in longer range 
correlations in the residual stresses. As an extreme ex- 
ample, one could consider the stress caused by segments 
of dislocations placed randomly in the solid with the con- 
straint that they form loops, i.e., that the dislocation 
density tensor is divergence free. This results in a = 1 
which is a strongly relevant perturbation yielding, in the 
linear approximation an unphysical roughness exponent 
of £h = 1 both for the front and the mode I crack surface. 
Generally, we have Oi = (3 — a)/2 for a < 3, thus very 



long range correlations with aw 1.4 would be needed to 
explain the experiments. 

On the basis of our results an explanation of the mea- 
sured roughness in terms of quasi-static motion of cracks 
appears unlikely. Only if very long range correlations in 
residual stresses existed would a quasi-static explanation 
be viable. However, even if these did occur, one would 
have replaced the problem of understanding the appar- 
ent universality of crack roughness exponents in a wide 
variety of materials with the problem of why residual 
stress correlations should be long range and universal. 
A more appealing alternative is that elastodynamics of 
the medium plays an essential role. It has been shown 
that, in the case of a crack front restricted to a plane, the 
sound waves emitted as it moves changes its behaviour 
both when it is moving at a finite velocity and near the 
threshold ju|. Such effects may also play a crucial role 
in increasing the roughness of the crack surface . 

We finally note that in thin plates, our quasi-static 
analysis yields a crack path with a roughness exponent 
of ( = 1/2 for tensile cracks and £ = 1 for tearing cracks; 
in these situations elastodynamic effects may be less im- 
portant, but buckling may play a role. 
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